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bJQi ABSTRACT 

' Context. The smoothed particle hydrodynamics (SPH) technique is a well-known numerical method that has been 

I applied to simulating the evolution of a wide variety of systems. Modern astrophysical applications of the method rely 

00 , on the Lagrangian formulation of fluid Euler equations, which is fully conservative. A different scheme, based on a 

■ matrix approach to the SPH equations is currently being used in computational fluid dynamics (CDF). An original 
matrix formulation of SPH based on an integral approach to the derivatives, called lADo, has been recently proposed 
and is fully conservative and well-suited to simulating astrophysical processes. 

Aims. The behavior of the lADo scheme is analyzed in connection with several astrophysical scenarios, and compared 
to the same simulations carried out with the standard SPH technique. 

Methods. The proposed hydrodynamic scheme is validated using a variety of numerical tests that cover important topics 
in astrophysics, such as the evolution of supernova remnants, the stability of self-gravitating bodies, and the coalescence 
of compact objects. 

^ ■ Results. The analysis of the hydrodynamical simulations of the above-mentioned astrophysical scenarios suggests that the 

( ■ SPH scheme built with the integral approach to the derivatives improves the results of the standard SPH technique. In 

I particular, there is a better development of hydrodynamic instabilities, a good description of self-gravitating structures 

, in equilibrium and a reasonable description of the process of coalescence of two white dwarfs. We also observed good 

conservations of energy and both linear and angular momenta that were generally better than those of standard SPH. 
In addition the new scheme is less susceptible to pairing instability. 
(N- Conclusions. We present a formalism based on a tensor approach to Euler SPH equations that we checked using a variety 

^ , of three-dimensional tests of astrophysical interest. This new scheme is more accurate because of the re-normalization 

(N. imposed on the interpolations, which is fully conservative and less prone to undergoing the pairing instability. The 

■ analysis of these test cases suggests that the method may improve the simulation of many astrophysical problems with 
\l ' only a moderate computational overload. 
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. 1. Introduction astrophysics could improve the simulations with the SPH 

technique with an affordable computational cost. In this 

. ^ . Multidimensional numerical hydrodynamics is one of the work, we focus on specific three-dimensional (3D) astro- 

K> . most powerful tools of modern astrophysics to comprehend physical applications of the scheme formulated in paper I. 

5_i ■ the cosmos machinery. Among them, the smoothed par- In particular, we choose examples from different fields of 

] tide hydrodynamics (SPH) is one of the most widely used astronomy to check the method and highlight its potential 

techniques because of its ability to describe the evolution of advantages over the standard SPH scheme. The suitability 

fluids with complicated geometries and a diversity of length of IADq for describing hydrodynamic instabilities found in 

sca les. Since it was formulated, m ore t han th i rty ye ars ago, paper I is confirmed by simulating the evolution of a su- 

by iGingold k Monaehanl (|l977t l and iLucvl (|1977D . it has pernova remnant (SNR). As the SNR evolves embedded in 

largely evolved incorporating, little by little, a plethora of an uniform background of particles with negligible gravity, 

methods that makes it competitive compared to grid-based there are no numerical troubles affecting the outer limits of 

methods of Eulerian type. Details of the modern mathemat- the system. The existence of boundaries becomes relevant 

ical formulation, as well as of the main features of the state- to the second test, which is devoted to describing the equi- 

of-art of the SPH te chnique, c an be found in the reviews librium features of polytropes with different indexes and 

bvlMonaghanl(|2005D .lR osswo g| (1200911 . Isl^ingell (l2010ll . and masses. For this problem, the interplay between pressure 



iPricd (|2012[ ) forces and gravity becomes crucial to ensure that the ob- 



iGarcia-Senz ^Tdl (|20T1 (henceforth p aper I) rece ntly Gained structures are compatible with the analytical mod- 
suggested that the use of matrix methods, [mti (fl999h . in els. We show that the tensor method leads to polytropes 
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where the central density and radius are slightly closer to 
the theoretical predictions than those obtained with the 
standard SPH technique. We also explore the ability of 
IADq to describe a very dynamical situation by simulat- 
ing the coalescence of two white dwarfs. In this 
catastrophic merging of the stars ensues after a few orbital 
periods. For this test, the tensor method gives results of, 
at least, comparable quality to those obtained using the 
standard SPH scheme, but displaying a more homogeneous 
mixing of the material of both stars. There are also some 
differences in the angular velocity distribution of the rem- 
nant. For a similar elapsed time, the matrix calculation does 
not lead to the complete rigid rotation of the core, which 
is, however, achieved in the simulation using the standard 
scheme. 

The text is organized as follows. In Sec. [H we review 
the mathematical formalism linked to IADq and discuss 
its most relevant features. In Sec. |3l we describe the three 
astrophysical tests aimed at validating the code and com- 
paring its performance to that of standard SPH. Section 0] 
is devoted to incorporating thermal conductive transport 
in the tensor scheme, and to check the resulting algorithm. 
The benchmarking of the code is done in Sec. [51 Finally, the 
main conclusions of our work, as well as some comments 
about the shortcomings of the developed scheme and fu- 
ture lines of improvement, are outlined in the last section, 
which is devoted to our conclusions. 



2. Main features of the lADn scheme 



where 



= X! ~ ^i,a){Xj.b - Xj,a)Wab{ha) ;i, j = 1,3 

b 

(5) 

and 



Ik,a = '^mb {xk,b - Xk,a)Wab{ha) = 1,3. (6) 
6 

It was shown in paper I that Eq. (|4]) leads to a formu- 
lation of the SPH Euler equations that is compatible with 
the variational principle. These equations are given by 



Pa = ^"ibWafc(|rb -r^\,ha), 



(7) 



(8) 



^) = ^ ^ni6(Wi,a-Ui,6) ( o^°„2 -^i.abiha) + Ai,ab 



dt 

where Ai^ab and -4 are 



pi 



(9) 



In paper I, it was shown that a conservative SPH scheme 
can be deduced from an integral approach to the deriva- 
tives. As starting point, we define the integral. 



Ai.ab{ha) = '^Cij,a{ha){Xj^b " Xj^a)Wab{ha) , (10) 



/(r)- / [/(r')-/(r)](r'-r)T4^(|r'-r|,/i)d, 
Jv 



,'3 



(1) 



where M^(|r' — r|, ft,) is a spherically symmetric interpolating 
function and h is called the smoothing length. The integral 
/(r) can be used to find the gradient of a function /(r) in a 
way similar to that used to evaluate the Laplace operator 
based on another integral expression in standard SPH. See 
iBrookshawl ()1985( ) and Monaghan (2005). The IADq inter- 
pretation of SPH is the consequence of approaching Eq. ^ 
with summations along with two reasonable simplifications 



/(rt) - /(ra) ~ Vfa • (rb - ra) , 



(2) 



where a and b refer to neighboring particles with masses 
nia and rUb, respectively, and 



•^iMbi^b) = ^ Cij^b{hb){Xj.b - Xj,a)Wab{hb) , (11) 



being Cij the coefficients of the inverse matrix defined in 
Eq. @ and d the dimension of the space. The magnitude 
r2a = (1 — dp/dh wLa dWab/dh) accounts for the gradi- 
ent of the smoothing length. As usual, liab gives the viscous 
pressure due to the artificial viscosity (AV) 



H„ 



— tor Yab Vab < U, ^-^2) 

otherwise. 



and the remaining symbols have their usual meaning. The 
coefficient fj,ab is 



V^/(rb)(rb-ra)I^(|rb-ra|,ft. 



(3) 



is the corresponding discrete expression for Eq. ([T]). Note 
that, because the kernel is spherically symmetric, the factor 
/(ra) does not appear in the expression above. 
The direct application of Eqs. (P),©, and ^ to calculating 
the gradient of density, leads to a matrix equation 



(4) 



dp/dxi ' 




'Til 


Tl2 


TlZ ' 


-1 




dp/dx2 




T21 




■^23 
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dp/ 0x3 . 
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. T3I 


T32 


■''33 . 







hab^ab ' '^ab n\ 

P'°-b = — , „ „i 72 ■ i^3j 

^ab + 0-01 Kb 

To compute the viscous acceleration, the arithmetic mean 
of A is taken to be 



A,afc = ^ [Ai^abiK) +A!Lab{hb)] 



(14) 



Therefore Eqs. ([7]), ^ and ^ summarize the basis of the 
IADq formalism. 

It is worth noting the strong similitudes between the 
above equations and those of standard SPH. The main dif- 
ference is that vector expressions have been changed to 
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tensor relationships that encode the renormalization of the 
many summations that appear in the SPH technique. Any 
expression of standard SPH can indeed be made compatible 
with IADq by taking the kernel derivative as 



(15) 



where the coefficients Ai,ab{ha) are defined by Eq. ([TT 
If the matrix coefficients in Eq. (U)) are calculated analyti- 
cally, the matrix T becomes diagonal. In this case, it can be 
shown that for Gaussian kernels the standard and IAD de- 
scriptions are totally equivalent. Hereafter, the scheme ob- 
tained from the diagonal form of T is referred to as vector- 
lADo. 



2.1. Harmonic kernels 

All the simulations presented in this work were car- 
ried out using the so-called harmonic kernels devised by 
ICabezon et al.l ([2008), to perform the SPH interpolations. 
Being relatively unknown kernels with very interesting fea- 
tures, their use warrants some explanation. They are de- 
fined as 



W^{v,h) = ^s^nc-[^v); 



0<v <2. 



(16) 



where sinc(u) = n is the index of the kernel, and 

Bn is the normalization constant. By defining smc(O) = 1, 
the function sine is extended to an analytical function with 
compact support on the real axis. Because of its connection 
with spectral analysis, the function sinc{u) is of special rel- 
evance to signal analysis, from where it borrows its name. 
The profile of the kernel for several values of the leading in- 
dex n is shown in Fig.[l] where it can be seen that the profile 
of W^{v,h) closely matches that of the cubic spline. The 
behavior of W;^ {v, h) is therefore very similar to that of the 
cubic spline, with the advantage that its second derivative is 
continuous and can be differentiated many times. Switching 
the index to n = 2 gives the {v, h) kernel, with a profile 
close to that of the truncated Gaussian kernel, as shown 
in Fig. [TJ The implementation of the W^{v,h) family of 
kernels adds more flexibility to SPH, as one can, for exam- 
ple, take a different index n to handle the artificial viscosity 
terms in either the momentum and energy equations or the 
heat conduction equation, without changing the number of 
neighbors of the particle. They are also useful for avoiding 
the pairing instability (see Sec. 12. 2p because whenever two 
particles approach each other too closely, the index of the 
kernel can be increased to block the development of the 
instability. 

The normalization consta nts Bn for a large nu mber of 
values of n were calculated in ICabezon et al.l ()2008[ ). where 
a fitting analytical formula for i?„ was also provided. To 
increase the computational speed, it is recommended that 
the value of smc and its derivative be stored in a table 
as a function of w, (0 < w < 2), and that a linear Taylor 
expansion be used to calculate the value of smc. This allows 
a fast computation of Eq. (jl6|) and its derivative once the 
index n has been chosen. 



2.2. Pairing instability 

The simulations oriented to benchmark the different SPH 
schemes reported in this manuscript as well as in paper I, 
suggest that there is an additional advantage of the matrix 
method. In general, calculations carried out using IADq are 
less affected by the pairing instability. For a given interpo- 
lating kernel with spherical symmetry, there is a fiducial 
distance to the center at which the first derivative of 
the kernel reaches its maximum absolute value. Within this 
critical radius, a pair of neighboring particles feels an in- 
creasingly weaker repulsive force, which can lead to the ar- 
tificia l clumping of the particles (but see lDehnen fc HossamI 
(2012) for a different explanation of the origin of this insta- 
bility). The exact location of rp depends on the number of 
neighbors and the peculiarities of the kernel. For either the 
cubic-spline kernel or the harmonic kernel with index n=3, 
its location is rg = while for more centrally condensed 
kernels it is at yet smaller radii. 

The robustness of matrix methods to cope with pairing 
instability can be understood by confronting the analytical 
estimation of the gradient of the kernel within the standard 
framework, to the numerical value obtained using Eq. (jlSp . 
The numerical derivative of Wab for the particle a located 
at the center of a 2D lattice was evaluated using two differ- 
ent particle settings and harmonic kernel indexes. Figure [5] 
shows the value of VaWab at different distances from the 
kernel origin. The left column stands for a bcc-type particle 
setting in a two-dimensional square lattice, while the col- 
umn on the right represents a quasi-uniform distribution of 
particles, obtained by randomly perturbing the bcc distri- 
bution with a maximum perturbation amplitude of 0.3 A, 
where A is the lattice spacing. The harmonic kernel index 
was set to rt = 3 (thus, reproducing the cubic-spline ker- 
nel) and to n = 5 to see the effect of making the profile 
sharper. It can be seen that increasing n (keeping h con- 
stant) always raises the maximum of the kernel derivative, 
as expected for a more centrally condensed kernel, shifting 
the abscissa where the maximum is achieved closer to zero. 
For n = 3 and n = 5, the maximum is taken at r/h = 2/3 
and r/h — 0.504, respectively. The effect of changing the 
numerical scheme also has an impact on the kernel deriva- 
tive. According to Fig. [51 the maximum of the gradient 
of the kernel always appears closer to the origin than in 
the standard scheme, regardless of the number of neigh- 
bors. Moreover, the absolute value of the maximum is larger 
for IADq. Together, this means that the tensor scheme is 
less affected by pairing instability than the standard SPH 
scheme. 



2.3. Free surface conditions 

Handling boundaries is often a difficult point of the SPH 
technique. In static gaseous configurations, such as stars 
or planets, the radius of the object is the result of the 
careful balance between gravity and pressure forces. The 
gradient in the pressure is, however, not as accurately es- 
timated by SPH near the surface as in the star's interior. 
Therefore, the radius of the configuration is a magnitude 
not as well-defined as others once the mechanical equilib- 
rium is achieved. 

Unfortunately, the direct application of IADq to free 
surface boundaries does not solve the problem. The reason 
is that the magnitudes r.y given by Eq. ([SJ, which serve to 
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normalize the derivatives, are very sensitive to the presence 
of boundaries. Close to the edge of the system, the gradi- 
ent of pressure is overestimated and the equilibrium of the 
body is reached at a larger radius than in the standard 
(STD) scheme. The impact of the free boundaries on the 
Tij coefficients is well-illustrated in Fig. [31 which depicts the 
profile of these magnitudes along the polytropic structure 
discussed in Sec. 13.21 The numerical estimation of can 
be compared to their analytical value calculated using 

—Bnh^ J v^sinc"'{—v) dv ;i — 1,3 , 

which for n = 3 gives rf^ = = T33 — 0.2916 h^. We see 
that, with the exception of the surface layers, the analytical 
and the numerical values agree with an accuracy better 
than a 5%. These coefficients are, however, very sensitive 
to boundaries and rapidly decay when the first particles 
belonging to the surface enter the summations. This led 
some authors (e.g. lOeer et 311120071 ) to propose a practical 
algorithm for choosing the adequate value of th as 

^,,^(ru (Eq.ED /3>/3o, .-^g^ 
" 1 T^^ otherwise, ^ ' 

where /3 = Tu/rf^i and /3o — 0.95. In Sec. 13.21 we discuss 
the impact of the different approaches to tu in describ- 
ing the structure of polytropes with different indexes. We 
chose several values for /3o ranging from /3o = (fully ten- 
sor lADo) to /3o = cx) (vector -IADn). The values /3o = 0-90 
and /3o = 0.97 suggested by iQger et all (|2007D were also 
checked. Our main conclusion is that the /3o = fully ten- 
sor IADq scheme provides the best description of the struc- 
ture of these polytropes, although the numerical noise is 
larger than for the STD scheme. Interestingly, the vector- 
like approach with /3o = 00 give results as good as the STD 
scheme with the same computational cost. Nevertheless, 
an additional advantage of the full tensor scheme is that 
it improves the description of the growth of hydrodynamic 
instabilities, as suggested in the tests described in Sec. [3] 
and paper I. On another note, we have not seen any partic- 
ular advantage of using a hybrid scheme with /3o — 0-95 to 
handle self-gravitating bodies, which has the negative side 
of increasing the numerical noise. Furthermore, a value of 
/3 < 1 is not necessarily linked to a free boundary because 
it may appear, for example, in the presence of mildly or 
strong shock waves. For these reasons, we conclude that 
the renormalized IADq scheme summarized by Eqs. ([5]) to 
(jlip should be preferentially used to carry out simulations 
of astrophysical systems. 

3. Astrophysical tests 

A basic check of the new scheme was described in paper 
I. For the most part, the technique was tested in two di- 
mensions with particles located in ordered lattices and, 
sometimes, using particles with different masses. For these 
specific tests, the IADq scheme showed a good behavior 
because it was able to handle the Ray leigh- Taylor and 
Kelvin-Helmholtz instabilities better than standard SPH. 



Simulation of supersonic events (Sedov blast wave and wall- 
heating shock) were of similar quality, if not slightly better, 
than those computed using STD smoothed particle hydro- 
dynamics. The total energy was always more well-conserved 
when the tensor method was used. 

The only 3D calculation in paper I was that used to sim- 
ulate the evolution towards stability of a Sun-like polytrope. 
In that simulation, gravity was calculated usin g a multipo- 
lar expansion of the force (jHernauist fc Kat3 ([198 9')). The 
main conclusion of the test was that, although good equi- 
librium configurations were achieved in both methods, in 
the tensor method more energy is stored as numerical noise 
for the same elapsed time. 

To check the code in realistic astrophysical scenarios, we 
chose three simulations related to different hydrodynamic 
processes. These include the growth of the RT instability in 
supernova remnants, the stability of polytropes of different 
masses and polytropic indexes, and the study of the coales- 
cence process of a pair of compact stars. The calculations 
were carried out by choosing n = 3 in the harmonic ker- 
nels given in Eq. 1161 Unless explicitly stated, the number 
of neighbors was kept constant to nf, = 100 in all simu- 
lations. Precise details of the initial particle setting and 
physics included are presented in each subsection of the 
corresponding test. 

3.1. Hydrodynamics of a supernova remnant 

The evolution of supernova remnants has to be studied in 
more than one dimension to capture the fine details of their 
structure. Several things may contribute to cause the evo- 
lution to deviate from the spherical symmetry: the lack of 
symmetry of the exploding object that gives rise to the 
SNR, the inhomogeneities in the ambient medium (AM), 
or th e interaction of the remnant with cosmic rays (|Wang| 
l201l[ ). One of the physical phenomenae that was among 
the first to be studied as a source of spatial irregularities 
in SNRs was the Rayleigh- Taylor (RT) instability. In SNR, 
the RT instability may often develop in the region between 
the forward and reverse shocks, induced by the deceleration 
of the supernova shell. Pioneering 3D simulations of the im- 
pact of the RT ins tability on the evoluti on of the remnant 
were conducted bv lChevalier et aD(jl992[ ). by taking a small 
slice of the full domain to enhance the resolution. More re- 
cently, 3D si mulations encompa ssing the whole SNR were 
developed by IVigh et aD (|201lt) . To reproduce the growth 
and structure of the instability, a large amount of com- 
putational cell s are needed. Tw o-dimensional simulations 
carried out bv iDwarkadasI (|2000f ) suggest that a minimum 
of 300x300 mesh zones (in 2D) have to be used to describe 
the growth of the RT fingers. Therefore, the ability to sim- 
ulate the growth of the RT instability in SNR is a good test 
for 3D hydrocodes. 

All multi-dimensional simulations performed so far 
have been carried out using grid-based codes, generally of 
Eulerian type. In spite of its wide use in astrophysics, SPH 
has scarcely been applied t o the study o f SNR. One ex- 
ception was the study of iGarda-Senz et all ()2012f ). who 
used an axisymmetric SPH code to study the imprint of 
the secondary star on the geometry of the remnant result- 
ing from a Type la supernova explosion. It is worth noting 
that SPH is able to describe the growth of the RT instabil- 
ities in 3D, using a smaher num ber of particles than those 
suggested in IDwarkadasI (|2000l) {N ~ 1.5 10^ particles in- 
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stead of 2 — 3 10^ grid cells) because there is no waste of 
computational resources in the lowest density regions of the 
simulated domain. 

To perform the SNR simulation, we built an exponen- 
tial profile for the e.jecta following the prescriptions given in 
iDwarkadai ()2000D . The profile was mapped to a 3D set of 
equal-mass particles and allowed to interact with an homo- 
geneous ambient medium. A perfect-gas equation of state 
(EOS), where 7 = 5/3 was assumed. 

The density of the ejecta was 

Pej{t)=Aexp(^-^^ t-\ (19) 
where the constants A and Ve are 

A = 7.67 10« g cm-3 s^ ( ^) ' E^,^ , (20) 
1 / M ■ \ 

Ve = 2.44 10*^ cm s'^ E^^ — ^ . (21) 

\Mch J 

We took A/ej — Mch and E^i = 1 in the expressions 
above. The time t in Eq. (fT9|) was set to io = 2 • 10^ s, and 
the velocity profile at to was assumed to be homologous, 
with v{r) — r/to. With these choices, the initial size of the 
ejecta was ~ 1 pc and the density at the outer edge matched 
that of the AM medium (assumed here to be homogeneous 
and have a value pam — 2 • 10"^'* g.crn"^). The ID den- 
sity profile of the ejecta was mapped onto a 3D sample of 
particles spread across a sphere of radius 1 pc, according 
to the radial density profile given by Eq. and where 

the angular location of the particles was chosen at random. 
This random setting led to a blurred density profile, which 
was finally driven to the spherical symmetry by means of a 
tangential relaxation of the model using SPH. To do this, 
we allowed the particles to move under the action of pres- 
sure forces, although their movement was constrained to en- 
sure that the distance to the center remained constant. The 
homogeneous AM was simply reproduced by sampling the 
particles in a cubic-centered, bcc grid with a size of 10 pc. 
Using Nej = 257, 776 particles and Nam = 1, 206, 576 par- 
ticles, the mass ratio of particles belonging to the AM to 
those in the ejecta is rngj /itiam — 2.3. To check that such 
a mass contrast has a negligible impact on the results, we 
recalculated model Bi of Table [1] using a number of parti- 
cles for the ambient medium N'j^-j^j — 2.3Nam- There were 
no significant changes in the evolution of the model. 

The evolution of the SNR was simulated for different 
initial conditions using both the STD and the IADq codes, 
with the goal of analyzing their performance. 

A single mode perturbation in the initial radial-velocity 
field was seeded at i = yr according to the expression 

Avr = (5 (1 - exp(-(r/i?ej)^)) cos(n 0) , (22) 

where R^j — 1 pc, 9 is the azimuthal angle, 9 = cos~^(z/r), 
and the parameter 6 is close to the maximum perturbation 
velocity we wish to impose. The wavenumber was set to 
n = 12 in all models. 

The inspection of the results summarized in Table [1] re- 
veals that energy is always better conserved, by a factor 
~ 2, when the tensor method is used. Linear and angu- 
lar momentum conservations are at least as good as in the 
standard formulation, independently of either the absence 



(models Ai and A2) or the presence (models Bi and B2) 
of the velocity perturbation. Figure |4] presents a color map 
of density in the YZ plane at t = 698 yr for models Bi 
and B2 of Table [1] calculated by assuming S = 500 km.s~^ 
in Eq. (P^ . The combined effect of the radial velocity per- 
turbation and the anisotropics of the bcc lattice leads to 
the growth of the RT instability. At t = 698 yr, the devel- 
opment of the instability is already appreciably larger in 
the tensor calculation than in the STD scheme. The differ- 
ences between both calculations become more pronounced 
at t = 951 yr, when the forward shock reaches the limits of 
the system. These results agree with the main conclusions 
of paper I, in terms of the development of the RT instability 
in 2D stratified fluids inside a homogeneous gravitational 
field. 

Model B3 in Table [T] and the lower panels of Fig. 2] re- 
fer to vector-IADo. As we can see, the evolution of the RT 
instability is similar to that of the STD scheme but total 
energy is better conserved and the density map seems to 
be slightly cleaner. Therefore, we conclude that the renor- 
malization of the derivatives, which characterizes the fully 
tensor IADq method, makes it more suitable for describ- 
ing the evolution of hydrodynamic instabilities. The im- 
portance of renormalization is highlighted in Fig. [5l which 
depicts the profile of the elements tu and Tij calculated 
with Eq. ([5]). As we can see, there is a large deviation of tu 
from its analytical estimation through Eq. (jl7p . Similarly, 
the off-diagonal matrix elements are not equal to zero. The 
re-normalization of the derivative then helps the growth 
of the instability, confirming the main conclusions given in 
paper I. 

As a final test, we built an developed model for the 
supern o va ejecta using the stretc hed-grid method (jHerantI 
(|1992D . iGarcfa-Senz et all (|1998D V These models are la- 
beled Ci and C2 in Table [T] This kind of initial models 
displays almost perfect spherically symmetric density pro- 
files, but the price to pay is a large deformation of the 
particle lattice. Although no perturbation was seeded, the 
larger anisotropics of the grid soon lead to the growth of 
the RT instability in the region between the forward and 
reverse shocks. The results of the simulations are shown 
in Fig. The upper panels depict the density profile af- 
ter t — 698 yr of evolution, whereas the bottom ones show 
the radial velocity profile at the same elapsed time. We see 
that outside the region where the RT instability develops 
the profiles calculated using IADq show less dispersion than 
those computed with the STD scheme. In the RT unstable 
layer located between the forward and reverse shocks, the 
dispersion is larger, as expected. In addition, there is a fac- 
tor of two enhancement in the energy conservation when 
the tensor method is used (last two rows in Table [I}. 

3.2. Polytropes of index n—3/2 and n=5/2 in equilibrium 

Gaseous self-gravitating structures can usually be roughly 
approached using a polytrope with the appropriate in- 
dex. Because of their simplicity and known theoreti- 
cal properties, polytropes are often used to benchmark 
multidimensional hydroc odes (jSteinmetz fc Miilleil 119931 : 
iPrice fc Monaghaiill2007l ). We explored the abilities of the 
tensor scheme in reproducing the equilibrium structure of 
polytropes of index n — 3/2 and n = 5/2, and confronted 
the results with the predictions of standard SPH models. 
A mass of M = O.6M0 and a radius of i? = 8 ■ 10* cm 
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were assigned to the n = 3/2 polytrope, so that it could 
represent a stable white dwarf with central density pc ~ 
3.3 • 10^ g.cm"^. For the n = 5/2 polytrope, we assigned it 
a mass of M = 1.35Mq and a radius of 2 • 10^ cm, which are 
both representative of a massive white dwarf with central 
density pc = 1-9 • 10^ g.cm~^. Note that in this last case the 
mass is close to the Chandrasekhar-mass limit for a poly- 
trope with n = 3. This second structure is then much more 
unstable than the 0.6Mq star. The initial models were built 
by randomly spreading 2 • 10'' particles in 3D, according to 
the mass density profile of the polytrope. In a second step, 
the models were relaxed in the non-radial direction to en- 
hance their spherical symmetry. Afterwards, the particles 
were allowed to evolve freely in the 3D domain. The EOS 
was that of a perfect gas P = p{'-f — l)u, with 7 = 5/3, 
7/5 for n — 3/2, 5/2 polytropes, respectively. The main 
features of both polytropes are summarized in Table [5J 

Figure [7] represents the evolution of the central density 
and stellar radius for models E and F of TablejH To smooth 
fluctuations, the central density and the radius were aver- 
aged over the closest particles to the center and the surface, 
respectively. The central density of the less massive struc- 
ture converges to the analytical value after ten oscillatory 
cycles. The radius of the polytrope follows a similar evo- 
lution, although the value obtained with STD is ~ 90% 
of the theoretical estimation owing to the small number of 
particles used in the simulation. We see that despite the 
lower convergence rate of the tensor calculation the final 
value of Pc and R are closer to the correct value than for 
the STD calculation. The relative error in the central den- 
sity is er{pc) < 2% for the IADq calculation and a little 
higher, er{pc) — 5% for STD. For the more massive poly- 
trope, the differences between the SPH schemes become 
more accentuated, as shown in the bottom panels of Fig. [T] 
The fluctuation in the central density is much larger than in 
the previous case and the convergence towards the analyt- 
ical value is less rapid. Nevertheless, the tensor calculation 
provides a more accurate estimation of the central density 
of the equilibrium configuration. According to Fig. [71 the 
relative errors in the central density at the last calculated 
models are er(pc) — 33% and er(Pc) — 42% for the tensor 
and standard methods, respectively. It is worth noting that 
a lower particle clumping was observed when the matrix 
scheme was used, which may be the origin of the small dif- 
ferences in the value of the central density calculated with 
both methods. 

Figure \8\ shows the evolution of both polytropes calcu- 
lated with vector-I ADq. The central density is much more 
closely reproduced than for the STD scheme but the radius 
of the configurations differs considerably from the actual 
value. Therefore, we conclude that the best combination of 
central density and surface radius is achieved by the full 
tensor IADq scheme. 

The leftmost panels of Fig. [HI depict the evolution of the 
ratio of the kinetic to internal energies, Ek/Eint, of both 
polytropes. The evolution of the kinetic energy is that of 
a damped system when standard SPH is used, while the 
IADq calculation is less dissipative. Therefore, the amount 
of kinetic energy stored in the configuration after several 
evolution cycles remains larger for the matrix calculation. 
All this points to the renormalization of the derivatives as 
the main agent behind the increase in numerical noise. To 
clarify this point, we ran two simulations, models E3 and 
F3 in Table [H using the vector approach to IADq. The 



evolution of the kinetic energy is shown in the rightmost 
panels of Fig. [9] As we can see, the damping of the systems 
is even faster that in the STD scheme and the final kinetic 
energy is lower. 

The level of energy conservation at the last calculated 
model is correlated with the amount of kinetic energy. 
According to the last column of Table [21 energy is most ac- 
curately conserved for the vector approach to IADq . For the 
n = 3/2 polytrope, there is a better conservation for STD 
than for tensor IADq but for the most unstable 71 = 5/2 
polytrope the level of conservation is similar. 

3.3. Merging of two white dwarfs 

Because of its mcsh-free features, a substantial amount of 
applications of SPH have been devoted to the study of the 
coalescence of stellar objects. In this section we simulated 
the merging of two twin polytropes of index n=3/2 and 
mass 0.6Mq with both the IADq and STD schemes. In each 
calculation, the structure was that of the last calculated 
model corresponding to models Ei and E2 in Table [51 The 
EOS used was that of a perfect gas with 7 — 5/3. Both 
stars were put in a circular rigid-rotation orbit in the plane 
XY with radius r°^j, = 1.5 Rs, where Rg = 8,000 km is 
the theoretical surface radius of the polytrope, from the 
center of mass of the system. To enforce the coalescence, a 
braking force proportional to the velocity Ffc^a — ~k{t) • v 
was imposed during one revolution period P, according to 
the prescription 

where t is the elapsed time and P = 29.3 s. 

This scenario roughly accounts for the merging of twin 
white dwarfs described by an EOS dominated by the elec- 
trons in the partially degenerated regime. 

The results of the simulations are summarized in 
Figs. [TUl to [131 Figures [TUl and [TT1 represent a density color 
map of both stars in the orbital plane as obtained with 
IADq and STD schemes, respectively. On the whole, the 
behavior is rather similar, although the coalescence evolves 
more slowly for the matrix method. We have to keep in 
mind, however, that the detai ls of the evolutio n depend on 
the precise initial conditions (|Dan et al.ll201lfl . The initial 
equilibrium models for both calculations are similar, but 
not identical, because of the differences introduced dur- 
ing the relaxation process. The model Ei of Table ^ ap- 
proached with the tensor method, has a larger radius and 
less clumping than model E2 calculated using the STD 
scheme. The details of the coalescence are also infiuenced 
by the errors in the estimation of gravity introduced by the 
tree-walk approach. These errors induce small differences 
that grow during the most dynamic phase of the merg- 
ing. Thus, it is not straightforward to discern which part 
of the differences comes from the particular SPH scheme 
used in the simulation. On the other hand, both meth- 
ods give the correct radius for Roche-Lobe (RL) overfiow 
once the center of mass separation, di2 , becomes similar to 
di2 — 17, 100 km. Once the RL is crossed, the mass transfer 
becomes catastrophic and the merging of both stars takes 
place in a small fraction of one period. 

Figure [T^l shows the evolution of the component of 
the angular momentum, orthogonal to the orbital plane. 
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For t < P, there is a constant loss of angular momentum 
due to the imposed external braking force. At times t > P, 
but before the dynamical phase of the merging {t/P ~ 4), 
the tensor scheme displays a rather good conservation of 
Lz, with small ripples induced by the multipolar approach 
to gravity. Conservation of in this phase is not so good 
for the standard formulation. On another note, during the 
most violent phase of the coalescence and further bouncing 
of the core, the evolution of Lz follows a slightly different 
trend in both schemes. The tensor method leads to a mono- 
tonic increase in Lz until a maximum is reached followed 
by a steady decline. In the STD calculation, the behavior 
is the opposite. Nevertheless, the evolution of the angular 
momentum after the merging is complicated by the forma- 
tion of an extended halo of particles, which is clearly seen in 
the last snapshot in Figs. [TU]and[TT] Although only a tiny 
amount of mass belongs to the halo, such mass produces 
a large torque on the core and any error in the estimation 
of gravity translates into a larger error in the angular mo- 
mentum. 

An important challenge for the hydrocodes is to ade- 
quately represent the mix of the advected material during 
the coalescence process. In our system, the initial condi- 
tions are fully symmetric. Therefore, one would expect that 
a few minutes after the merging the core is homogeneously 
composed of the material of both stars. Figure [T3] depicts 
the approximate distribution of the gas belonging to each 
star one minute after the catastrophic stage of the coales- 
cence process. As we can see, the material of both white 
dwarfs is much more thoroughly mixed in the IADq cal- 
culation than the standard one. The merged core is made 
of successive thin, onion-like layers. In the calculation with 
the standard scheme, the onion-like structure is less pro- 
nounced, especially in the inner 10 km. We stress that the 
recipe to handle the artificial viscosity is exactly the same 
in both calculations. The only difference is the treatment 
of the gradient of the kernel, which seems ultimately to be 
the responsible for the larger amount of mixing obtained in 
the IADq calculation. 

Figure [14] shows the angular velocity profile of the par- 
ticles in the orbital plane as a function of their mass coor- 
dinate. As we can see, the standard calculation leads to an 
almost perfect rigid rotation for Air < 0.8Mq, followed by 
a Keplerian velocity distribution beyond tha t point. This 
behav ior agrees with the results obtained by iRaskin et al.l 
(|2012f) for a similar scenario. At distances Mr > IMq, 
the velocity profile obtained with IADq is also Keplerian, 
but below that mass, rigid rotation is never attained in the 
matrix calculation. In the standard calculation, the high 
amount of viscosity, which prevents the mixing of the core, 
now couples the different layers of the fiuid so that the 
system rapidly approaches rigid rotation. Nevertheless, the 
time delay for core synchronization in nature is a function 
of the real physical viscosity. In the hydrodynamic models, 
the synchronization time is probably artificially shortened 
by the much larger numerical viscosity introduced by the 
codes. The impact of artificial viscosity in the fin al structure 
of the merged object was extensively discussed in lDan et al.l 
(l20Tll) . 

4. Adding thermal conduction to IADq 

It is not difficult to incorporate physics into the proposed 
matrix formulation of SPH by simply using Eq. ([T5|) to 



calculate the gradient of the kernel whenever necessary. As 
an example, we work out an expression that can be used 
to calculate the conductive heat transport and apply it to 
simulate the evolution of a thermal wave. To do this, it is 
enough to cons ider the SPH thermal conduction equation 
(lSDringelll20Tol) 



du, 
~dt 



a 



nib (Ka + Kb){Tb - To) 



6=1 



Pa Pb 



Vab ■ VaVFab , (24) 



' ab 



where k is the thermal conductivity. If we estimate the gra- 
dient of the kernel using Eq. ([T5|) . it leads to, 



dug 
dt 



nib {Kg + Kb){Tb - Tg) 



6=1 



Pa Pb 



^ ^ (^^^,0 '^i^b) ^i,ab ; 



1=1 



(25) 

where d is the dimension of the space and the tilde symbol 
means the arithmetic average of the magnitude, Eq. (ITi)) . 

We applied Eq. (^5)) to the classical test of simulat- 
ing the propagation of a thermal wave in an homogeneous 
medium with p = 1 g.cm~^ and compared the results 
with those obtained using Eq. . During the calculation, 
we keep the particles at rest so that the energy equation, 
Eq. dni), reduces to the heat transport equation given by 
the expressions above. An initial point-like discontinuity in 
energy was seeded at the center of the box. To avoid nu- 
merical problems, the discontinuity was smoothed using a 
sharp Gaussian with characteristic width of a few times 
the smoothing length parameter. To achieve this we made 
use of the analytical e xpression for a spherically symmetric 
thermal- wave front in iLandau fc Lifshit3 (jl959tl 



u{r,t) = 



A 



{Ana t) 2 



exp 



4a t 



uo , 



(26) 



where a = K/{cy p) is the thermal diffusivity and Ci, the 
specific heat capacity. The parameters in Eq. (j26|) were set 
to A = 10^ erg cm^ g~^, uq = 10'^ erg g^^,c„ = 9.1 • 
106 erg g-i K-i, and k = 3.9 • 10^ erg s"^ cm"! K'^ 

To set the initial model, we put 47, 707 identical parti- 
cles at the nodes of a 3D regular lattice. The initial energy 
profile was that given by Eq. (j26l) for t = 0.04 s. The evolu- 
tion of the thermal wave for different numbers of neighbors 
is depicted in Fig. (TSl As we can see, the tensor scheme is 
in closer agreement with the analytical solution when the 
number of neighbors is small, Uf, ~ 25. Nevertheless, above 
Ub ~ 50 neighbors both schemes lead to similar results. 
To investigate the influence of the initial setting of parti- 
cles we carried out a second calculation, this time using a 
pseudo-ordered lattice. The new particle distribution was 
obtained after randomly perturbing the 3D ordered lattice 
with a maximum perturbation amplitude 0.3A, where A 
is the lattice spacing. In this new calculation the mass of 
the particles was conveniently arranged to keep the density 
constant, at p = 1 g cm~^, throughout the system. The 
maximum mass contrast between neighboring particles was 
lower than a factor of two. The results of this second cal- 
culation are summarized in Fig. 1161 As we can see, the 
calculation with IADq always reproduces more closely the 
evolution of the heat wave, especially for a small number of 
neighbors but also even when a moderate, Ub ~ 50, amount 
of neighbors is chosen. 
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5. Computational issues 

In this section we conducted a benchmark test to assess the 
impact of the calculation of the IADq terms and the mo- 
mentum and energy equations using Eqs. ((Sl llll) . To do this 
we generated a 3D distribution of particles using a quasi- 
random Sobol distribution. Using the stretching technique, 
we then displaced the particles to follow a spherical density 
profile for th e central core of a pre-collapse 15 Mq star (see 
iHeeer et all (2()fl5)). The hydrocode is OMP-parallelized, 
and includes a Barnes-Hut octal-tree solver to calculating 
the gravitational force up to the quadrupolar terrn , and the 
Lattimer-Swesty EOS (jLattimer fc Swestvl (llQQlD l. As we 
are only interested in the impact of the IADq calculation 
on the overall wall-clock time, we set the velocity of the 
particles to zero at each time-step (frozen structure). 

Table [3] summarizes the outcome of the calculations. We 
varied the number of particles and the number of threads 
to gain insight into both the overhead and the scaling of 
the IADq calculation. Table [3] shows, from left to right, the 
number of particles, the number of threads, the tolerance 
parameter 9, the wall-clock time per iteration for the stan- 
dard calculation (i™^), the wall-clock time spent by the 
code in the calculation of the standard momentum and en- 
ergy equations (t^^^^), the wall-clock time per iteration for 
the IADq calculation {I^aDq ) > wall-clock time spent by 
the code in the calculation of the IADq coefficients (tf^l;^), 
and the wall-clock time spent by the code in the calculation 
of the IADq momentum and energy equations {tf^AD^)- The 
last column shows the IADq overhead as a percentage, with 
respect to the standard simulation, and is calculated as 



IADq overhead 



iCAL I 4-MOM 



t 



MOM 
STD 



t 



TOT 
STD 



XlOO. (27) 



All these times were averaged over 1000-2000 iterations for 
each calculation, which was always performed on the same 
12-core AMD machine. 

The parameter 9 was used to control the criterion for 
choosing between opening a node in the tree (and using 
the particle-to-particle interaction to calculate gravity) and 
not opening the node (and using its global contribution via 
the multipolar expansion). The smaller the value of 9, the 
slower the tree-walk, because more nodes are opened and 
the calculation gets closer to N'^ scaling. Typical values 
of 6*, those actually used in calculations, are between 0.5 
and 0.7. Here, we selected two values (0.3 and 0.6) to as- 
sess the relevance of the lADg-related calculations to the 
gravity evaluation, which is typically the most important 
bottleneck of the code. 

From these results, it can be seen that the inclusion 
of the IADq scheme in single-threaded (serial) calculations 
contributes with an overhead of approximately 18%, inde- 
pendently of the number of particles used, which means 
that, as expected, it has a linear dependence on them. This 
can be seen in Fig. [T7] (left), where it is clear that the in- 
clusion of IADq in the calculation, does not vary the slope 
of the curve. In the case of the standard calculation, the 
relationship between the wall-clock time and the number 
of particles is log{t) ~ logiN)^-"^"'' , while for the lADo 
calculation it is log{t) ~ log{N)^-°^^^ . Even more, Fig. [TT] 
(right) shows that the overhead ratio of the lADo-related 
calculations to their standard analogous tends to saturate 



at a factor slightly lower than 2. Noting that the momen- 
tum and energy calculations in the standard version of SPH 
takes around 20% of the time, this result is consistent with 
the IADq overhead. 

These are desirable properties as they mean that this 
section of the code will have at least the same behavior in 
terms of parallelization as the rest of the code, and that the 
overhead is limited independently of the number of parti- 
cles. This can also be seen when we compare the single- 
threaded results with the multi-threaded ones. The wall- 
clock time for the IADq sections in the 12-thread calculation 
is 9-10 times shorter than the time for the same sections in 
the serial work, which points to a good strong scaling and 
a considerable reduction in the IADq overhead to 9% of the 
overall wall-clock time. 

The overhead reduction due to parallelization of course 
depends on how well the code is parallelized and the amount 
of it that remains serial. Nevertheless, the results that refer 
specifically to the IADq sections show relevant information 
on the low impact that IADq may have on the overall cal- 
culation when these sections are parallelized. 

Finally, it can also be seen from the last two rows of 
Table [3] that the relevance of the IADq sections is reduced 
to ~ 14% when gravity becomes more dominant. This can 
be taken as a lower limit of the IADq overhead in single- 
threaded calculations, as the parameter 9 = 0.3 is rather 
extreme and is rarely used below that value. 

6. Conclusions 

We have checked the behavior of a novel SPH scheme where 
gradients are calculated using an integral approach. The 
main features of the technique, called IADq, were described 
in detail in paper I and summarized in Sec. [5] of this pa- 
per. The main virtue of the approach relies in that the re- 
normalization of the derivatives appears naturally, without 
any degradation of the conservative properties that charac- 
terize the SPH technique. Another relevant feature is that 
the basic mathematical formalism of IADq looks very simi- 
lar to that of a standard SPH technique, making its imple- 
mentation quite straightforward. As commented in paper I, 
matrix methods based on a similar, although not identical, 
formulation have b een us ed in CDF for the past decade (see 
for instance Diltij (jigggf )). but have never been previously 
applied to astrophysics. 

Three test cases of considerable astrophysical interest 
were selected to validate IADq, as well as to detect its 
virtues and weaknesses. The performance of the method 
in describing the growth of the RT instability in a super- 
nova remnant was analyzed in Sec. 13.11 with the conclu- 
sion that IADq provides a healthier development of the RT 
fingers. The method's success in handling hydrodynamical 
instabilities relies on the re-normalization imposed on the 
derivatives, thus confirming the results obtained in paper I 
using toy models. 

The second test was addressed especially to the appli- 
cations of the method to describing stellar objects (ap- 
proached as polytropes with different indexes and known 
analytical properties). One important question here relies 
in the treatment of the outer boundary of the o bject, which 
is often a controversial point in matrix methods (I Peer et al.l 
2007,). We have entirely explored this issue using several 
recipes to handle the surface of the polytropes. These in- 
clude the use of tensor and vector formulations of IADq , as 
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well as hybrid schemes that use the full matrix expressions 
in the interior but changes to vector-IADo near the surface, 
according to Eq. (|18l) . The best results were obtained when 
the full tensor IADq scheme was used, especially for the 
most unstable polytrope with index n=5/2. Nevertheless, 
the amount of numerical noise stored as residual kinetic 
energy at equilibrium was larger for tensor-IADo than for 
STD. 

The ability of the tensor method to handle the co- 
alescence and further merging of stellar-like objects was 
checked in Sec. 13.31 Even though for this particular test 
the simulations did not show any clear advantage of the 
matrix method, it gave a good depiction of the coalescence 
process. Despite both schemes being conservative by con- 
struction, complete preservation of angular momentum is 
impeded by the use of the hierarchical cluster method in 
calculating gravity. Nevertheless, there are indications that 
the angular momentum-component orthogonal to the or- 
bital plane behaves better in IADq than in the standard 
scheme (Fig. [T^ . The prompt product after the coalescence 
is a core surrounded by an extended diluted halo of parti- 
cles moving at Keplerian velocity. The calculations show 
that the tensor method leads to a more homogeneous mix- 
ing of the material of both stars in the core region. As the 
numerical recipe for implementing AV in both calculations 
is exactly the same, we conclude that the different behavior 
is caused by the differences in the algorithm used to com- 
pute the kernel derivative. Therefore, it seems that IADq 
is able to provide a more effective hydrodynamical mixing 
than the standard scheme, but still avoid the penetration 
of fluids in strong shocks, as suggested in paper I. These 
differences also affect the distribution of angular velocity 
in the core of the remnant, which in the STD calculation 
soon reaches rigid rotation but in the matrix one does not, 
for a similar elapsed time. A potential weakness of the ten- 
sor method in handling these kinds of configurations comes 
from the computation and further inversion of matrix T 
calculated with Eq. ([5]). In the case of the isolation of a par- 
ticle, as it could be for those belonging to the most diluted 
region of the domain, matrix 7" becomes singular leading 
to the complete breakdown of the simulation. In this sense, 
the matrix scheme is less robust than the standard one. 
Nevertheless, current SPH schemes usually include an al- 
gorithm to ensure that the number of neighbors of a given 
particle remain constant during the simulation. This algo- 
rithm is necessary to reliably compute the magnitude VL in 
Eqs. (|S]) and (O and, when working properly, to avoid the 
singularity problems linked to matrix 7". 

Another interesting features of the matrix method are 
its ability to avoid the pairing instability and the better 
handling of thermal conduction. The first one leads to less 
particle clustering, thus improving the quality of the inter- 
polations; while in the second case, the results of Sec. |3]sug- 
gest that diffusion-like equations can be handled by IADq 
in a better way than in the standard framework. 

We also conducted a benchmark test to evaluate the 
impact of the inclusion of the IADq formalism on the wall- 
clock time for a nominal calculation. According to Table [3] 
and Fig. [T71 the computational overload introduced in a 
serial calculation by the re-normalization of the derivatives 
is low (~ 20%), being practically independent of the to- 
tal number of particles. The scaling of the code with the 
number of particles remains virtually untouched, and the 



lADp-related sections of the code show a good behavior in 
front of parallelization. 

Therefore, the analysis of the astrophysical tests dis- 
cussed in this paper support the main conclusions of paper 
I, where the basic implementation of the integral approach 
to the derivatives was discussed and checked. All this sug- 
gests that the use of the re-normalized, fully conservative 
IADq approach to the SPH equations may improve, in gen- 
eral, the quality of the simulations in astrophysics with a 
little computational overload penalty. 

The smaller amount of viscous dissipation shown by 
IADq, compared to STD for the same AV formalism, sug- 
gests that it could be of interest not only when handling 
hydrodynamic instabilities but also when simulating tur- 
bulence. Turbulence is at the heart of many astrophysi- 
cal problems, b eing of especial relevan ce to understanding 
star formation ()Federrath et al.l ()201Cl[ )). In this respect, a 
full comparison of 3D astrophysical turbulence calculated 
with a variety of algorithms, both SPH and grid codes, 
was provided by 'Kitsionas et alj (|2009[ ) , iPrice fc Federratbl 
(2013), andlKritsuk et al. (20l!| , with the result that both 
families of codes give similar results for an equivalent 
number of resolution elements in each direction in space. 
Nevertheless, these experiments also show that, owing to 
the higher dissipation, the scaling range of SPH codes is 
slightly sh orter than that of grid -based codes, as demon- 
strated bv iKitsionas et all ()2009[ ). Therefore, it may be of 
great interest to check the ability of the proposed lADg 
scheme to handle astrophysical turbulence in the near fu- 
ture. 
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Fig. 1. Profile of several kernels in one dimension. The 
Gaussian kernel is truncated air — 2h. The Wn{H) belongs 
to the harmonic family of kernels described by Eq. (1161) for 
n — 2, 3, and 6, respectively. 
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Fig. 2. Profile of the gradient of the kernel calculated an- 
alytically in the standard SPH way (referred as analytical 
in the figures) as well as numerically, through Eq. ([T5|) for 
two particle settings: hcc and pseudo-random, kernel in- 
dexes n = 3, 5, and different number of neighbors. 
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Fig. 3. Profile of coefficients Ta/h? and T.ij/h'^,{j ^ i) in 
a polytrope of index n — 5/2 fitted with 20,000 particles 
(model Fi of Table [U . The distance to the center is nor- 
malized to the radius of the polytrope, Rg. The analytical 
value of Tii/h? is also given (horizontal dashed line corre- 
sponding to /3 = 1). Basically, all the star interior has a 
value /3 > 0.95, while beyond a radius Rs ~ 2hs (where hg 
is the smoothing length close to the surface) the value of /3 
rapidly declines (see Sec. 12.31 for the definition of /?). 
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Fig. 4. Density color map depicting the growth of the RT 
instabihty in a SNR for models Bi (upper panels), B2 (mid- 
dle panels), and B3 (lower panels) of Table [T] at times 
t — 698 yr (left colmrms) and t = 951 yr (right colmxins). 
The size of the box in all panels is 10 pc in each direction. 
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Fig. 5. Profile of the coefficients Ta/h? and Tij/h? ^ {j ^ i) 
for model Bi of Table [1] The analytical value of Tu/h"^ is 
also given (horizontal dashed line). 
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Fig. 6. Density and velocity profiles of the SNR at elapsed 
time t = 698 yr for the stretched-grid models of the su- 
pernova ejecta (models Ci and C2 in Table [1]). Outside the 
RT-unstable region, the spherical symmetry is better pre- 
served when the tensor method is used. 
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Fig. 7. Evolution of central density and surface radius (nor- 
malized to the theoretical values) for polytropes with in- 
dexes n — 3/2 (upper panels, corresponding to models Ei 
and E2 in Table [5]) and n — 5/2 (lower panels, correspond- 
ing to models Fi and F2 in Table[2]). Continuum and dashed 
lines are for the tensor-lADo and standard schemes, respec- 
tively. 
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Fig. 8. Evolution of central density and surface radius (nor- 
malized to the theoretical values) for polytropes with in- 
dexes n = 3/2 (upper panels, corresponding to models E2 
and E3 in Table [2]) and n = 5/2 (lower panels, correspond- 
ing to models F2 and F3). Continuous and dashed lines are 
for vector-IADo and standard schemes, respectively. 
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Fig. 9. Evolution of the ratio of kinetic to internal energies. 
Left column is for polytropes with indexes n = 3/2 (upper 
panel) and n = 5/2 (lower panel) calculated using tensor- 
IADq (continuum line) and STD (dashed line). Same for 
pictures on the right column, where, however, vector-IADo 
was used. 
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Fig. 10. Density color map of the coalescence process of 
two twin polytropes of index n — 3/2 at times t = 0.31P, 
t = 2.8F, t = 4.3F, and t = 7.7P, {P = 29.3 s) calculated 
using IADq. 




Fig. 11. Density color map of the coalescence process of 
two twin polytropes of index n = 3/2 at times t = 0.23P, 
t = 2.15P, t = 3.3P, and t = 7.4P (P = 29.3 s) calculated 
using the standard SPH scheme. 
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Fig. 12. Left: Evolution of the distance between the center 
of mass of both polytropes. Right: Evolution of the orbital 
angular momentum, L^-, normalized to its initial value. An 
artificial braking force was acting before t/P = 1 




Fig. 13. Mixing of material of both stars in the core as 
calculated by IADq at t = 227 s (upper picture) and with 
STD a.t t = 216 s (bottom picture). Blue and red colors 
refer to the gas belonging to each component of the orig- 
inal binary system. Mixed regions display a superposition 
of both colors. 
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Fig. 14. Angular velocity profiles at t ~ 220 s as function 
of mass coordinate, calculated using IADq (red) and STD 
(blue) , respectively. The Keplerian values for both calcula- 
tions are given for reference. 



Fig. 16. Evolution of a thermal wave propagating through 
a pseudo-ordered distribution of particles at different times. 
The profile of internal energy is shown for both SPH 
schemes and two different numbers of neighbors. 
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Fig. 15. Evolution of a thermal wave propagating through 
an ordered lattice of particles at different times. The profile 
of internal energy is shown for both SPH schemes and two 
different numbers of neighbors. 
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Fig. 17. Benchmarking of the impact of IADq on the com- 
putational time. Left: Total wall-clock time per time-step 
in front of the number of particles. Right: Overload ratio 
of the IADq related calculations (matrix coefficients + mo- 
mentum and energy equations) to the standard momentum 
and energy calculation. 
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Model Scheme Particles Ejecta Perturbation S | AE\/Eo \Arcm\/R \L\/iJ2i \Li\) 







(xlO**) 




(km.s ^) 








Ai 


lADo 


1.46 


rnd 





6.5 10"^ 


8 IQ-'^ 




A2 


STD 


1.46 


rnd 





1.3 10"^ 


8.2 IQ-^ 


4 10"* 


Bi 


lADo 


1.46 


rnd 


500 


6.5 10"^ 


7.5 IQ-^ 


5.5 10"^ 


B2 


STD 


1.46 


rnd 


500 


1.3 10~^ 


8.4 10"^ 


4 10-* 


B3 


lADo (vector) 


1.46 


rnd 


500 


7 10"^ 


7.5 IQ-^ 


3 10"* 


Ci 


lADo 


1.46 


Stretch 





6.5 10"^ 


9 10"^^ 


4 10-12 


C2 


STD 


1.46 


Stretch 





1.3 10"^ 


7 10"" 


9 10"^ 



Table 1. Comparison between several magnitudes at i = 698 yr of evolution of the SNR computed using IADq and STD 
schemes. The deviation of the center of mass is normalized to R = 2.5 pc. \Li\ = {\Lxi\ + \Lyi\ + \Lzi\). 



Model scheme Mass (Mq) index n pe (g-cm ") (cm) Apc/pc ARs/Rs Ekin/Etot \ AE\/Eo 



El 


lADo 


0.6 


3/2 


3.3 10** 


8 lO** 


1 10"^ 


0.02 


6 10-* 


2 10-^ 


E2 


STD 


0.6 


3/2 


3.3 10^ 


8 10* 


6 10-2 


0.11 


3 10-^ 


4 10-* 


E3 


lADo (vector) 


0.6 


3/2 


3.3 10'^ 


8 10*^ 


3 10-2 


0.14 


4 lO"-""' 


2 10-* 


Fi 


lADo 


1.35 


5/2 


1.9 10^ 


2 10** 


0.33 


0.04 


8 10-* 


2 10-^ 


F2 


STD 


1.35 


5/2 


1.9 10^ 


2 10* 


0.42 


0.11 


2 10-* 


2 10-^ 


F3 


lADo (vector) 


1.35 


5/2 


1.9 10^ 


2 10* 


5 10-^ 


0.24 


5 10-^ 


3 10-* 



Table 2. Main fcatm-cs of polytropcs with polytropic indexes n = 3/2 and n = 5/2, respectively. The analytical value 
of the central density pc and surface radius Rs arc given in columns 4 and 5. The relative errors between the numerical 
and analytical estimations axe provided in columns 7 and 8. The kinetic energy stored as numerical noise (normalized to 
the internal energy) and the level of energy conservation at the last calculated model are given in columns 9 and 10. 



Particles 


Threads 


9 


'■STD 


J.MUM 
''STD 


'■IADq 


'■lADn 


J.MUM 
'■lADn 


lADo overhead 


20,000 


1 


0.6 


3.84 


0.85 


4.51 


0.55 


0.99 


17.97% 


50,000 


1 


0.6 


11.27 


2.45 


13.32 


1.40 


3.05 


17.75% 


150,000 


1 


0.6 


40.29 


8.12 


47.18 


4.80 


10.26 


17.22% 


500,000 


1 


0.6 


148.87 


28.24 


175.40 


18.11 


36.42 


17.66% 


1,500,000 


1 


0.6 


474.68 


88.50 


558.69 


56.33 


116.56 


17.78% 


5,000,000 


1 


0.6 


1750.95 


341.29 


2082.26 


198.11 


473.78 


18.88% 


500,000 


12 


0.6 


27.27 


3.73 


30.37 


1.73 


4.33 


8.50% 


5,000,000 


12 


0.6 


311.51 


38.08 


347.76 


18.23 


48.70 


9.26% 


50,000 


1 


0.3 


14.28 


2.47 


16.46 


1.42 


3.09 


14.29% 


150,000 


1 


0.3 


48.25 


8.16 


55.10 


4.81 


10.22 


14.24% 



Table 3. Average wall-clock time for the different benchmark tests and sections of the code relevant to the IADq 
calculation. 6 is the parameter for opening the nodes in the tree when calculating gravity. TOT, MOM, and CAL stand 
for total time per iteration, momentum and energy equation calculations, and IADq terms calculation, respectively. All 
times are expressed in seconds. 
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